#include "NeutrinoTools.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TDirectory.h"

#include "style.h"

TString getTag(const TString fin)
{
  TString tag;
  if(fin.Contains("Neut")){
    tag="Neut";
  }
  else  if(fin.Contains("Genie")){
    tag="Genie";
  }
  else  if(fin.Contains("NuWro")){
    tag="NuWro";
  }
  else{
    printf("getTag wrong file name %s\n", fin.Data()); exit(1);
  }

  return tag;
}

void plotPID(const TString fin)
{
  TFile::Open(fin);
  const TString tag = getTag(fin); 

  TTree * tt=(TTree*) gDirectory->Get("tree");

  TCanvas * cc=new TCanvas;
  gStyle->SetTitleW(0.95);
  gStyle->SetStatX(0.35);
  gStyle->SetStatY(0.8);
  gStyle->SetOptStat("nemrou");
  gStyle->SetOptTitle(1);

  gPad->SetLogy();

  TH1D * hh= new TH1D("hh","",4, -0.5, 3.5);
  hh->SetMarkerSize(3);
  hh->SetMaximum(100000);
  hh->SetMinimum(0.9);

  const TString var = "(fProtonSimPDG==2212)*2+(fMuonSimPDG==13)";
  const TString cut0 = "fProtonNRec>0 && fMuonNRec>0 && fNeutrinoType==14 ";

  hh->SetXTitle(var);
  hh->SetYTitle("counts");

  hh->SetName(tag+"CCQE");    hh->SetTitle(cut0+" &&  (fNeutMode==1||fNeutMode==2)"); 
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sPID.png", hh->GetName()));

  hh->SetName(tag+"bk");      hh->SetTitle(cut0+" && !(fNeutMode==1||fNeutMode==2)"); 
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sPID.png", hh->GetName()));
}

void plotTPC()
{
  TCanvas * cc=new TCanvas;

  TTree * tt=(TTree*) gDirectory->Get("default");

  TH2D *h2=new TH2D("h2","",100,200,3000,200,0,1000);
  NeutrinoTools::BinLog(h2->GetXaxis());

  TH2D *h3=(TH2D*)h2->Clone("h3");

  tt->Draw("selp_tpc_dedx_raw:selp_mom>>h2",  "accum_level[][1]==11","colz");
  tt->Draw("selmu_tpc_dedx_raw:selmu_mom>>h3","accum_level[][1]==11","colz");
  h2->Draw("colz");
  h3->Draw("same");

  cc->Print("outplot/tpc.png");

}

void plotStat(const TString fin)
{

  TFile::Open(fin);
  const TString tag = getTag(fin); 

  TTree * tt=(TTree*) gDirectory->Get("tree");

  TCanvas * cc=new TCanvas;
  gStyle->SetTitleW(0.95);
  gStyle->SetStatX(0.35);
  gStyle->SetStatY(0.8);
  gStyle->SetOptStat("nemrou");
  gStyle->SetOptTitle(1);

  gPad->SetLogy();

  TH1D * hh= new TH1D("hh","",14,11.5, 15.5);
  hh->SetMarkerSize(3);
  hh->SetMaximum(100000);
  hh->SetMinimum(0.9);

  //-------------------- Type 0 
  TString var = "fProtonNRec*8+fMuonNRec*4+fProtonNSim*2+fMuonNSim";   
  hh->SetXTitle(var);
  hh->SetYTitle("counts");

  hh->SetName(tag+"Type0all");       hh->SetTitle("fEventType == 0");
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sstat.png", hh->GetName()));

  hh->SetName(tag+"Type0nonQBS");    hh->SetTitle("fEventType == 0 && fNeutrinoType==14 && fNeutMode==1"); 
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sstat.png", hh->GetName()));

  hh->SetName(tag+"Type0QBS");       hh->SetTitle("fEventType == 0 && fNeutrinoType==14 && fNeutMode==2"); 
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sstat.png", hh->GetName()));

  hh->SetName(tag+"Type0NumuBK");    hh->SetTitle("fEventType == 0 && fNeutrinoType==14 && (fNeutMode!=1 && fNeutMode!=2)"); 
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sstat.png", hh->GetName()));

  hh->SetName(tag+"Type0nonNumuBK"); hh->SetTitle("fEventType == 0 && fNeutrinoType!=14"); 
  tt->Draw(Form("%s>> %s", var.Data(), hh->GetName()), hh->GetTitle(),"text"); cc->Print(Form("outplot/%sstat.png", hh->GetName()));
}

void plotFM()
{
  TCanvas * c0= new TCanvas;

  TTree * tt=(TTree*) gDirectory->Get("tree");
  tt->Draw("fNeutronToy->P()","fMuonNSim>0");
  c0->Print("outplot/FMp.png");

  tt->Draw("cos(fNeutronToy->Theta())","fMuonNSim>0");
  c0->Print("outplot/FMtheta.png");

  tt->Draw("fNeutronToy->Phi()","fMuonNSim>0");
  c0->Print("outplot/FMphi.png");

}

Bool_t IsRatio(const Int_t dopt)
{
  return dopt<=0;
}

void ptBroadenVarDraw(TTree * tmptree, TString tag, const TString var, const TString vname, const TString cut, const Int_t dopt, 
                      TCanvas * c0, TCanvas * cc, TLegend *ll, const Int_t nx, const Double_t xmin, const Double_t xmax, const Double_t non0start,  const TString title, const Int_t d1, const Int_t col1, const Int_t msty, const Int_t lsty, const Double_t hmax)
{
  //===========================================================================================
  //controls
  //===========================================================================================

  TString sRef;
  TString sType;
  TString detCut; 
  if(tag.Contains("TOY")){
    sRef="Toy";
  }
  else if(tag.Contains("SIM")){
    sRef="Sim";
    //sType=" vtx";
  }
  else if(tag.Contains("DET")){
    sRef="Det";
    //sType=" det";
    if(tag.Contains("PID")){
      detCut=" && fProtonSimPDG==2212 && fMuonSimPDG==13 ";
    }
  }
  else{
    printf("ptBroadenVarDraw wrong tag %s\n", tag.Data()); exit(1);
  }

  const Bool_t kREC = tag.Contains("REC");

  printf("test sRef %s kREC %d %s\n", sRef.Data(), kREC, tag.Data());
  
  const Bool_t kratio = IsRatio(dopt);
  if(kratio&&!kREC){
    printf("ptBroadenVarDraw wrong kratio %d kREC %d %s\n", kratio, kREC, tag.Data());exit(1);
  }

  tag += Form("v%sk%d",var.Data(), dopt);
  tag.ReplaceAll("-","_");

  const TString ctit=Form(";%s;%s", vname.Data(),kratio?"scaled to sim.":"#frac{d#it{N}}{#it{N}d#it{x}}");

  //===========================================================================================
  //drawing
  //===========================================================================================

  c0->cd();

  TH1D * href = new TH1D(tag+"ref","", nx, xmin, xmax); style::ResetStyle(href);   
  TString refdvar;
  if(var=="Mag"||var=="Theta"||var=="Phi"){
    refdvar=Form("fProton%sDeltaPt->%s()>> %s", sRef.Data(), var.Data(), href->GetName());
  }
  else if(var.Contains("Kt")){
    if(var.Contains("Proton"))
      refdvar=Form("fProton%sKt>> %s", sRef.Data(), href->GetName());
    else if(var.Contains("Muon"))
      refdvar=Form("fMuon%sKt>> %s", sRef.Data(), href->GetName());
  }

  const TString refdcut=cut+detCut+" && fProtonNSim>0 && fMuonNSim>0";
  if(non0start>0){
    NeutrinoTools::BinLog(href->GetXaxis(),non0start);
  }
  const Int_t refent = tmptree->Draw(refdvar,refdcut);
  printf("check refdvar | %50s |, refdcut | %s |, refent %d\n", refdvar.Data(), refdcut.Data(), refent);
  href->SetTitle(ctit);

  href->SetLineColor(col1);
  href->SetLineWidth(lsty==kSolid?2:4);
  href->SetLineStyle(lsty);

  ll->SetHeader("CCQE a.o. (sim)");

  TH1D * hcmp = 0x0;
  if(kREC){
    hcmp = new TH1D(tag+"cmp","", nx, xmin, xmax); style::ResetStyle(hcmp);
    TString recdvar;
    if(var=="Mag"||var=="Theta"||var=="Phi"){
      recdvar=Form("fProtonRecDeltaPt->%s()>> %s", var.Data(), hcmp->GetName());
    }
    else if(var.Contains("Kt")){
      if(var.Contains("Proton"))
        recdvar=Form("fProtonRecKt>> %s", hcmp->GetName());
      else if(var.Contains("Muon"))
        recdvar=Form("fMuonRecKt>> %s", hcmp->GetName());
    }

    const TString recdcut=cut+detCut+" && fProtonNRec>0 && fMuonNRec>0";
    if(non0start>0){
      NeutrinoTools::BinLog(hcmp->GetXaxis(),non0start);
    }
    const Int_t recent = tmptree->Draw(recdvar,recdcut);
    printf("check recdvar | %50s |, recdcut | %s |, recent %d\n", recdvar.Data(), recdcut.Data(), recent);
    hcmp->SetTitle(ctit);

    hcmp->SetMarkerStyle(msty);
    hcmp->SetMarkerSize(msty==29?3:2);
    hcmp->SetMarkerColor(col1);
    hcmp->SetLineColor(col1);
    hcmp->SetLineWidth(href->GetLineWidth());

    ll->SetHeader("CCQE a.o. (#bullet rec/-sim)");
  }

  cc->cd();

  if(kratio){    //kREC=1
    NeutrinoTools::ScaleToRef(hcmp, href);
    NeutrinoTools::ScaleToRef(href, href);
    
    //hcmp->SetMaximum(var=="Phi"?1.5:2);
    //hcmp->SetMinimum(var=="Phi"?0.5:0);
    
    hcmp->SetMaximum(1.5);
    hcmp->SetMinimum(0.5);
    
    href->SetFillColor(kBlack);
    href->SetFillStyle(0);
  }
  else{
    TH1D * htmp =0x0;
    htmp = NeutrinoTools::ToPDF(href); delete href; href = htmp;
    href->SetMinimum(0);
    href->SetMaximum(hmax);
    
    if(kREC){
      htmp = NeutrinoTools::ToPDF(hcmp); delete hcmp; hcmp = htmp;
      hcmp->SetMinimum(0);
      hcmp->SetMaximum(hmax);
    }
  }

  if(kREC){
    hcmp->Draw(d1? "same e":"e");
    href->Draw(kratio?"same E3":"same hist");

    ll->AddEntry(hcmp, title+sType, "p");
  }
  else{
    //href->SetFillColor(col1);
    //href->SetFillStyle(0);
    href->Draw(d1?"same hist e":"hist e");
    ll->AddEntry(href, title+sType, "l");
  }

  //if(var!="Theta"){
  if(var=="Phi"){
    gPad->SetLogx();
    if(dopt>=6 && dopt<=11){
      href->SetMinimum(1e-1);
      gPad->SetLogy();
    }
  }

  if(kratio){
    TF1 * f1 = new TF1("f1","1",0,10);
    f1->SetLineWidth(1);
    f1->SetLineColor(kBlack);
    f1->SetLineStyle(kDashed);
    f1->Draw("same");

    TF1 * f2 = new TF1("f2","1.05",0,10);
    f2->SetLineWidth(1);
    f2->SetLineColor(kRed);
    f2->SetLineStyle(kDashed);
    f2->Draw("same");

    TF1 * f3 = new TF1("f3","0.95",0,10);
    f3->SetLineWidth(1);
    f3->SetLineColor(kRed);
    f3->SetLineStyle(kDashed);
    f3->Draw("same");
  }
}

void ptBroadenSubDraw(TTree * tmptree, const Int_t dopt, TCanvas * ccs[], TLegend * lls[], const TString tag, const TString title, const TString cut, const Int_t d1, const Int_t col1, const Int_t msty, const Int_t lsty, const Double_t ymax[])
{
  ptBroadenVarDraw(tmptree, tag, "Mag",   "#delta p_{T} (GeV/c)",           cut, dopt, ccs[0], ccs[1], lls[1], 15, 0, 0.75,    0.05, title, d1, col1, msty, lsty, ymax[1]);
  ptBroadenVarDraw(tmptree, tag, "Phi",   "#delta#phi_{T} (rad.)",         cut, dopt, ccs[0], ccs[2], lls[2], 10, 0.01, 3.15, 0.05, title, d1, col1, msty, lsty, ymax[2]);
  ptBroadenVarDraw(tmptree, tag, "Theta", "#delta#alpha_{T} (rad.)",       cut, dopt, ccs[0], ccs[3], lls[3], 20, 0, 3.15,    -999, title, d1, col1, msty, lsty, ymax[3]);
  ptBroadenVarDraw(tmptree, tag, "ProtonKt", "k_{T}^{p} (GeV/c)", cut, dopt, ccs[0], ccs[4], lls[4], 15, 0, 0.75,    0.05, title, d1, col1, msty, lsty, ymax[4]);
  ptBroadenVarDraw(tmptree, tag, "MuonKt", "k_{T}^{#mu} (GeV/c)",             cut, dopt, ccs[0], ccs[5], lls[5], 15, 0, 0.75,    0.05, title, d1, col1, msty, lsty, ymax[5]);

  ccs[0]->cd();
  
  TH2D * h2dsim = new TH2D(tag+"h2dsim","",20, 0,3.142,20, 0,3.142);
  tmptree->Draw(Form("fProtonSimDeltaPt->Theta():fProtonSimDeltaPt->Phi()>>%s", h2dsim->GetName()), cut);

  TH2D * h2drec = new TH2D(tag+"h2drec","",20, 0,3.142,20, 0,3.142);
  tmptree->Draw(Form("fProtonRecDeltaPt->Theta():fProtonRecDeltaPt->Phi()>>%s", h2drec->GetName()), cut);

  //---  //---

  ccs[6]->cd();

  style::ResetStyle(h2drec); style::ResetStyle(h2dsim);
  h2dsim->SetTitle(title+" line: simulation, point: reconstruction;#delta#phi_{T} (rad.);#delta#alpha_{T} (rad.)");
  h2dsim->SetLineColor(col1);
  //h2drec->SetMarkerStyle(msty);
  //h2drec->SetMarkerSize(1);
  h2drec->SetMinimum(0);
  //h2drec->Draw(d1? "same cont1":"cont1");
  h2dsim->Draw(d1? "same box" : "box");
}

TString getQACut(const TString tag)
{
  if(tag=="dphiT"){
    const TString cmuon   = "fMuonRecFlightPath->Pt()>1.5   && 0.25<(fMuonRecChi2/fMuonRecNDOF)     && (fMuonRecChi2/fMuonRecNDOF)<4  && ";
    const TString cproton = "fProtonRecFlightPath->Pt()>0.5 && 0.25<(fProtonRecChi2/fProtonRecNDOF) && (fProtonRecChi2/fProtonRecNDOF)<4  && ";
    return cmuon + cproton;
  }
  else if(tag=="dalphaT"){
    const TString cmuon   = "fMuonRecNhits[0]>60  && ";
    const TString cproton = "fProtonRecNhits[0]>60  && ";
    return cmuon + cproton;
  }
  else{
    printf("getQACut bad tag %s\n", tag.Data()); exit(1);
  }
}

void ptBroadening(const Int_t dopt, const TString fneut, const TString fgenie, const TString fnuwro)
{
  const Int_t kratio = IsRatio(dopt);

  style::SetGlobalStyle();

  TTree * tneut=0x0;
  TTree * tgenie=0x0;
  TTree * tnuwro=0x0;

  if(fneut!="xxx"){
    TFile * finneut = new TFile(fneut);
    tneut = (TTree*)finneut->Get("tree");
  }
  if(fgenie!="xxx"){
    TFile * fingenie = new TFile(fgenie);
    tgenie = (TTree*)fingenie->Get("tree"); 
  }
  if(fnuwro!="xxx"){
    TFile * finnuwro = new TFile(fnuwro);
    tnuwro = (TTree*)finnuwro->Get("tree"); 
  }
  
  style::fgkYTitleOffset = 1.4;
  style::fgkXTitleOffset = 1.4;

  const Double_t lm=0.15;
  const Double_t rm=0.02;
  const Double_t tm=0.03;//12
  const Double_t bm=0.17;

  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gStyle->SetTitleW(0.95);
  TCanvas * c0=new TCanvas("c0","",800,600); 

  const Double_t lx = 0;//kratio?-0.48:0;
  const Double_t ly = kratio?-0.5:0;

  TCanvas * c1=new TCanvas("c1","",800,600);
  style::PadSetup(c1);
  c1->SetLeftMargin(lm);
  c1->SetRightMargin(rm);
  c1->SetTopMargin(tm);
  c1->SetBottomMargin(bm);
  TLegend * l1=new TLegend(0.65+lx, 0.67+ly,0.96+lx,0.94+ly);
  style::ResetStyle(l1,0.12);
  //l1->SetBorderSize(1);
  //l1->SetHeader("CCQE rec (-- sim)");

  TCanvas * c2=new TCanvas("c2","",800,600);
  style::PadSetup(c2);
  c2->SetLeftMargin(lm);
  c2->SetRightMargin(rm);
  c2->SetTopMargin(tm);
  c2->SetBottomMargin(bm);
  TLegend * l2=new TLegend(0.65+lx, 0.67+ly,0.96+lx,0.94+ly);
  style::ResetStyle(l2,0.12);
  //l2->SetBorderSize(1);
  //l2->SetHeader("CCQE rec (-- sim)");

  TCanvas * c3=new TCanvas("c3","",800,600);
  style::PadSetup(c3);
  c3->SetLeftMargin(lm);
  c3->SetRightMargin(rm);
  c3->SetTopMargin(tm);
  c3->SetBottomMargin(bm);
  TLegend * l3=new TLegend(0.32, 0.67+ly,0.63,0.94+ly);
  style::ResetStyle(l3,0.12);
  //l3->SetBorderSize(1);
  //l3->SetHeader("CCQE rec (-- sim)");

  TCanvas * c4=new TCanvas("c4","",800,600);
  style::PadSetup(c4);
  c4->SetLeftMargin(lm);
  c4->SetRightMargin(rm);
  c4->SetTopMargin(tm);
  c4->SetBottomMargin(bm);
  TLegend * l4=new TLegend(0.32, 0.67+ly,0.63,0.94+ly);
  style::ResetStyle(l4,0.12);
  //l4->SetBorderSize(1);
  //l4->SetHeader("CCQE rec (-- sim)");

  TCanvas * c5=new TCanvas("c5","",800,600);
  style::PadSetup(c5);
  c5->SetLeftMargin(lm);
  c5->SetRightMargin(rm);
  c5->SetTopMargin(tm);
  c5->SetBottomMargin(bm);
  TLegend * l5=new TLegend(0.32, 0.67+ly,0.63,0.94+ly);
  style::ResetStyle(l5,0.12);
  //l5->SetBorderSize(1);
  //l5->SetHeader("CCQE rec (-- sim)");

  TCanvas * c6=new TCanvas("c6","",800,600);
  style::PadSetup(c6);
  c6->SetLeftMargin(lm);
  c6->SetRightMargin(rm);
  c6->SetTopMargin(tm);
  c6->SetBottomMargin(bm);

  TCanvas * ccs[]={c0,  c1, c2, c3, c4, c5, c6};
  TLegend * lls[]={0x0, l1, l2, l3, l4, l5};
  //===========================================================================
  //===========================================================================

  //don't care about numu or non-numu, the contamination is at per mil level
  //non-QBS and QBS both have nearly all Nsim>0, absolutely all Nrec>0
  //for bk, has to apply PID cut so that Nsim and Nrec is comparable
  const TString qacut = "1 &&";// getQACut("dalphaT"); //"1 &&";// 

  const TString sigcut = qacut + "1";

  const TString bkcut = qacut + "(fNeutMode!=1 && fNeutMode!=2)";

  TTree * tmptree = 0x0;
  TString gen;

  if(kratio){
    if(dopt==-1){
      tmptree = tneut;
      gen="Neut";
    }
    else if(dopt==-2){
      tmptree = tgenie;
      gen="Genie";
    }
    else if(dopt==0){
      tmptree = tnuwro;
      gen="NuWro";
    }
    //dummy
    const Double_t ymax[]={100, 5.9, 7, 0.78, 10, 10};
    //ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC",gen+" signal", sigcut+" && (fNeutMode==1||fNeutMode==2)",0, kBlack,  20, kDashed, ymax);
    //ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC",gen, sigcut,0, kBlack,  20, kDashed, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC",gen+" true #mu p", sigcut+" && (fProtonSimPDG == 2212 && fMuonSimPDG == 13)",0, kBlack,  20, kDashed, ymax);
  }
  if(dopt==1){
    const Double_t ymax[]={100, 4.9, 3.5, 0.9, 10, 10};//1.28};
    tmptree = tneut;
    gen="Neut";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC1",gen+" non-QBS",      sigcut+" && fNeutMode==1",0, kRed,    21, kSolid, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC2",gen+" QBS",          sigcut+" && fNeutMode==2",1, kBlue,   20, kDashed, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC3PID",gen+" bk (true #mu p)",           bkcut,                    1, kGreen+3,29, kDotted, ymax);
  }
  else if(dopt==2){
    const Double_t ymax[]={100,5.5, 8.5, 0.91, 10, 10};
    tmptree = tgenie;
    gen="Genie";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC1",gen+" non-QBS",      sigcut+" && fNeutMode==1",0, kRed, 21, kSolid, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC2PID",gen+" bk (true #mu p)",           bkcut,                    1, kBlue,20, kDashed, ymax);
  }
  else if(dopt==3){
    const Double_t ymax[]={100,6.4, 4.9, 0.82, 10, 10};
    tmptree = tnuwro;
    gen="NuWro";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC1",gen+" non-QBS",       sigcut+" && fNeutMode==1",0, kRed,    21, kSolid, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC2",gen+" QBS",           sigcut+" && fNeutMode==2",1, kBlue,   20, kDashed, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIMREC3PID",gen+" bk (true #mu p)",            bkcut,                    1, kGreen+3,29, kDotted, ymax);
  }
  else if(dopt==4){
    //cmp all CCQE1
    const TString tmpcut = sigcut+" && fNeutMode==1";
    const Double_t ymax[]={100,6.6, 8.1, 0.65, 10, 10};
    
    const TString sRef="SIM";

    tmptree = tneut;
    gen="Neut";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+sRef,gen+" non-QBS", tmpcut,0, kGreen+3, 21, kSolid, ymax);
    
    tmptree = tgenie;
    gen="Genie";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+sRef,gen+" non-QBS", tmpcut,1, kRed, 29, kDotted, ymax);
    
    tmptree = tnuwro;
    gen="NuWro";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+sRef,gen+" non-QBS", tmpcut,1, kBlue, 20, kDashed, ymax);
  }
  else if(dopt==5){
    //cmp all CCQE2
    const TString tmpcut = sigcut+" && fNeutMode==2";
    const Double_t ymax[]={100,3.6, 1.7,  0.75, 10, 10};
    
    const TString sRef="SIM";

    tmptree = tneut;
    gen="Neut";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+sRef,gen+" QBS", tmpcut,0, kRed, 21, kSolid, ymax);
    
    tmptree = tnuwro;
    gen="NuWro";
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+sRef,gen+" QBS", tmpcut,1, kBlue, 20, kDashed, ymax);
  }
  else if(dopt==6 || dopt==7 || dopt==8){
    const Double_t ymax[]={100,11.5, 13.5, 0.82, 10, 10};
    if(dopt==6){
      tmptree = tneut;
      gen="Neut";
    }
    else if(dopt==7){
      tmptree = tgenie;
      gen="Genie";
    }
    else {
      tmptree = tnuwro;
      gen="NuWro";
    }
    
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"TOY",    "toy F-M only",       sigcut+" && fNeutMode==1",                       0, kRed,29, kSolid, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIM1",   gen+" non-QBS",       sigcut+" && fNeutMode==1",                       1, kGreen+3, 21, kDotted, ymax);
    if(gen!="Genie"){
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIM2",   gen+" QBS",           sigcut+" && fNeutMode==2",                       1, kBlue,20, kDashed, ymax);
    }
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIM3PID",   gen+" bk (true #mu p)",            bkcut,                                       1, kMagenta,23, kDashDotted, ymax);
  }
  else if(dopt==9 || dopt==10 || dopt==11){
    printf("obsolete!\n"); exit(1);

    const Double_t ymax[]={100,6.3, 8.1, 1.28, 10, 10};
    if(dopt==9){
      tmptree = tneut;
      gen="Neut";
    }
    else if(dopt==10){
      tmptree = tgenie;
      gen="Genie";
    }
    else{
      tmptree = tnuwro;
      gen="NuWro";
    }
    
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIM1",   gen+" non-QBS",   sigcut+" && fNeutMode==1",0, kRed, 21, kDashed, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"DET1",   gen+" non-QBS",   sigcut+" && fNeutMode==1",1, kRed, 22, kSolid, ymax);
    if(gen!="Genie"){
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"SIM2",   gen+" QBS",       sigcut+" && fNeutMode==2",1, kBlue,20, kDashed, ymax);
    ptBroadenSubDraw(tmptree, dopt,ccs,lls,gen+"DET2",   gen+" QBS",       sigcut+" && fNeutMode==2",1, kBlue,20, kSolid, ymax);
    }
  }
  
  c1->cd(); l1->Draw(); c1->Print(Form("outplot/ptBroadeningMag_kratio%d_opt%d.png",kratio, dopt));
  c2->cd(); l2->Draw(); c2->Print(Form("outplot/ptBroadeningPhi_kratio%d_opt%d.png",kratio, dopt));
  c3->cd(); l3->Draw(); c3->Print(Form("outplot/ptBroadeningTheta_kratio%d_opt%d.png",kratio, dopt));
  c4->cd(); l3->Draw(); c4->Print(Form("outplot/ptBroadeningPTT_kratio%d_opt%d.png",kratio, dopt));
  c5->cd(); l3->Draw(); c5->Print(Form("outplot/ptBroadeningMuonKT_kratio%d_opt%d.png",kratio, dopt));
  c6->cd(); c6->Print(Form("outplot/ptBroadening2d_kratio%d_opt%d.png",kratio, dopt));
}

//==========================================================================================================================
//==========================================================================================================================
